x <- 1:400
y <- sin(x / 10) * exp(x * -0.01)
plot(x,y)

philly <- readOGR("Philly3")
## OGR data source with driver: ESRI Shapefile 
## Source: "Philly3", layer: "Philly3"
## with 381 features
## It has 12 fields
names(philly)
##  [1] "OBJECTID"  "AREA"      "PERIMETER" "TRACTID"   "FIPSSTCO" 
##  [6] "TRT2000"   "STFID"     "N_HOMIC"   "TOTALPO"   "mdHHnc_"  
## [11] "HOMIC_R"   "PCT_COL"
plot(philly)

spplot(philly)

# RColorBrewer for colors
library(RColorBrewer)
display.brewer.all(type="seq")

pal <- brewer.pal(5, 'OrRd')
typeof(pal) #character
## [1] "character"
philly$HOMIC_R
##   [1]  426.74253   22.70148  122.45898   24.52784   22.70148   29.73536
##   [7]  122.45898   65.25711  122.45898   65.25711  122.45898  122.45898
##  [13]  406.12309   29.73536  122.45898   65.25711  122.45898   19.49698
##  [19]   42.58944   65.25711  595.23810    0.00000   22.41147   34.19973
##  [25]    0.00000  118.97680    0.00000    0.00000    0.00000    0.00000
##  [31]  126.18297  595.23810  141.37606   72.93946   22.41147   46.88233
##  [37]    0.00000  328.46715    0.00000   34.19973  118.97680    0.00000
##  [43]  190.80328   22.41147    0.00000    0.00000    0.00000    0.00000
##  [49]  190.80328   66.22517  328.46715   66.22517   46.88233   59.48840
##  [55]   19.49698   72.93946  190.80328  141.37606  595.23810    0.00000
##  [61]  259.74026  141.37606   66.22517    0.00000  269.98650  259.74026
##  [67]    0.00000   66.22517   66.22517  193.67334  259.74026   72.93946
##  [73]  259.74026   48.94763   95.58402  259.74026  299.53917  595.23810
##  [79]   66.22517  109.52903   40.76641    0.00000  259.74026  211.41649
##  [85]   40.76641   86.13264  193.67334   88.82309  259.74026   95.58402
##  [91]  595.23810   88.82309    0.00000  126.28541   86.13264  595.23810
##  [97]  259.74026  180.73845   95.58402  180.73845  211.41649   88.82309
## [103]   95.58402  245.48092  109.52903  595.23810   72.96607  111.11111
## [109]   86.13264   88.82309   40.76641  279.13468   70.12623  372.12449
## [115]  175.39464  100.83549  180.73845  134.37850  109.52903   45.46143
## [121]  465.11628  121.40132   86.13264   82.06138    0.00000  294.77821
## [127]  279.13468  100.83549  217.83526   45.46143  180.73845    0.00000
## [133]   45.46143  877.54288  877.54288  217.83526  253.36500  279.13468
## [139]  279.13468  502.18694  502.18694  180.73845    0.00000    0.00000
## [145]   45.46143  330.50047   45.46143    0.00000  285.19575    0.00000
## [151]  260.75619  306.69620  638.29787  285.19575  330.50047  877.54288
## [157]  253.67834  638.29787   45.46143  918.03279  502.18694    0.00000
## [163]  350.20694  877.54288  226.78594  285.19575  260.75619  877.54288
## [169]    0.00000 1343.87352  285.19575  226.78594  350.20694  663.50711
## [175]  136.96285  638.29787  328.52450  404.70225  244.89796  638.29787
## [181]   23.18572  136.96285  350.20694  663.50711  226.78594  656.25707
## [187]  117.78563  988.21741  244.89796  244.89796  663.50711  656.25707
## [193]  226.78594  288.73917  112.88805  426.49272  663.50711  456.10034
## [199]  240.32684   57.58157  456.10034  155.35992  567.06114  240.32684
## [205]  567.06114  112.88805  663.50711   53.03633  112.88805  663.50711
## [211]  663.50711  537.41215  537.41215  988.21741  755.12406  689.51931
## [217]  293.82958  302.43729  726.16538  726.16538   38.66976  293.82958
## [223] 1080.19640  539.32584 1067.86083  223.71365  237.60494  210.18216
## [229]   38.66976  726.16538  539.32584   54.79452  539.32584   53.03633
## [235]   38.66976  538.07272  203.09723  781.44694  764.08787  293.82958
## [241]  449.16292  210.18216  781.44694  163.20821  479.23323  781.44694
## [247]   54.79452  465.11628  347.03684  516.94428  293.82958  553.86320
## [253]  764.08787  137.42556  479.23323  210.18216   28.54696  600.24010
## [259]  347.03684  465.11628  465.11628  764.08787  764.08787  449.16292
## [265]  421.11173  236.62177  449.16292  764.08787  348.04909  347.03684
## [271]  764.08787  553.86320  137.42556   68.44627  280.11204  449.16292
## [277]   68.44627  301.84123  280.11204  245.26980  553.86320  240.21963
## [283]  354.05193  347.03684  312.06657  613.49693  613.49693  301.84123
## [289]  220.88987  663.34992  663.34992  220.88987  312.06657  287.52156
## [295]  312.06657  470.58824  392.43667  159.13431  380.16791  470.58824
## [301]  177.17034  177.17034  250.11034   27.70851    0.00000    0.00000
## [307]   38.41721   38.41721   27.70851  143.64376  143.64376    0.00000
## [313]  143.64376  177.17034    0.00000  104.42047  143.64376  250.11034
## [319]  345.78147  147.71049  354.47026  147.71049   40.46945  175.47280
## [325]  104.42047  354.47026  104.42047   74.57122   40.46945  343.81139
## [331]   74.57122  345.78147  175.47280  354.47026   74.57122  382.40918
## [337]  563.44377  354.47026   74.57122   74.57122  594.87179  880.00000
## [343]  115.63367  115.63367  608.64273  363.14117  354.47026  880.00000
## [349]   40.01601  117.39845  363.14117  240.26910  354.47026  240.26910
## [355]  880.00000   78.35455   54.47612  880.00000   40.01601   76.96254
## [361]  240.26910   54.47612  240.26910    0.00000  160.10247  124.22360
## [367]  360.03600   54.47612    0.00000   54.47612   54.47612    0.00000
## [373]    0.00000  304.41400    0.00000  304.41400    0.00000    0.00000
## [379]  304.41400  567.06114  567.06114
spplot(philly,"HOMIC_R", col.regions = pal, cuts = 4)  #applying colors and cuts

#===============================================================
# Denver datasets
denver1 <- readOGR("county_boundary_lines")
## OGR data source with driver: ESRI Shapefile 
## Source: "county_boundary_lines", layer: "county_boundary_lines"
## with 936 features
## It has 13 fields
plot(denver1)

# spplot(denver1)

denver2 <- readOGR("county_boundary")
## OGR data source with driver: ESRI Shapefile 
## Source: "county_boundary", layer: "county_boundary"
## with 12 features
## It has 1 fields
plot(denver2)

spplot(denver2)

# denver3 <- readOGR("tl_2017_us_county")
# plot(denver3)
# spplot(denver3)

denver4 <- readOGR("statisticalneighborhoods")
## OGR data source with driver: ESRI Shapefile 
## Source: "statisticalneighborhoods", layer: "statistical_neighborhoods"
## with 78 features
## It has 4 fields
plot(denver4)

spplot(denver4)

names(denver4)#"NBHD_ID"   "NBHD_NAME" "TYPOLOGY"  "NOTES"   
## [1] "NBHD_ID"   "NBHD_NAME" "TYPOLOGY"  "NOTES"
denver4$NBHD_NAME
##  [1] Chaffee Park                 Sunnyside                   
##  [3] Highland                     Globeville                  
##  [5] Jefferson Park               Sun Valley                  
##  [7] Valverde                     Athmar Park                 
##  [9] Windsor                      Northeast Park Hill         
## [11] Elyria Swansea               Wellshire                   
## [13] University                   Rosedale                    
## [15] Cheesman Park                Hilltop                     
## [17] Montclair                    Hale                        
## [19] North Park Hill              South Park Hill             
## [21] University Park              Platt Park                  
## [23] College View - South Platte  Overland                    
## [25] Ruby Hill                    Kennedy                     
## [27] Hampden                      Baker                       
## [29] Fort Logan                   Bear Valley                 
## [31] Harvey Park South            Southmoor Park              
## [33] Hampden South                Indian Creek                
## [35] Goldsmith                    Virginia Village            
## [37] Gateway - Green Valley Ranch DIA                         
## [39] University Hills             Harvey Park                 
## [41] Mar Lee                      Westwood                    
## [43] East Colfax                  Auraria                     
## [45] Cory - Merrill               Belcaro                     
## [47] Washington Park              Washington Park West        
## [49] Speer                        Cherry Creek                
## [51] Country Club                 Congress Park               
## [53] City Park                    Clayton                     
## [55] Skyland                      Cole                        
## [57] Marston                      Washington Virginia Vale    
## [59] Barnum                       Barnum West                 
## [61] Villa Park                   West Colfax                 
## [63] West Highland                Sloan Lake                  
## [65] Berkeley                     Regis                       
## [67] Lincoln Park                 City Park West              
## [69] Whittier                     Capitol Hill                
## [71] North Capitol Hill           Civic Center                
## [73] CBD                          Union Station               
## [75] Five Points                  Stapleton                   
## [77] Montbello                    Lowry Field                 
## 78 Levels: Athmar Park Auraria Baker Barnum Barnum West ... Windsor
# denver4$TYPOLOGY #NA
spplot(denver4, "NBHD_NAME", col.regions = pal, cuts = 78)

d5pal <- brewer.pal(8, 'YlOrRd')
denver5 <- readOGR("animal_care_and_control_division_boundaries")
## OGR data source with driver: ESRI Shapefile 
## Source: "animal_care_and_control_division_boundaries", layer: "animal_care_and_control_division_boundaries"
## with 8 features
## It has 4 fields
plot(denver5)

# spplot(denver5)
names(denver5)#"DISTRICT_I" "INSPECTOR_" "INSPECTOR1" "INSPECTO_1"  
## [1] "DISTRICT_I" "INSPECTOR_" "INSPECTOR1" "INSPECTO_1"
denver5$DISTRICT_I
## [1] F26 F27 F29 F32 F31 F33 F30 F28
## Levels: F26 F27 F28 F29 F30 F31 F32 F33
# denver5$INSPECTO_1 #INSPECTO_1
spplot(denver5, "DISTRICT_I", col.regions = d5pal, cuts = 78)

library(classInt)
breaks_qt <- classIntervals(philly$HOMIC_R, n = 5) # 5 colors
str(breaks_qt)
## List of 2
##  $ var : num [1:381] 426.7 22.7 122.5 24.5 22.7 ...
##  $ brks: num [1:6] 0 53 134 261 471 ...
##  - attr(*, "style")= chr "quantile"
##  - attr(*, "nobs")= int 156
##  - attr(*, "call")= language classIntervals(var = philly$HOMIC_R, n = 5)
##  - attr(*, "intervalClosure")= chr "left"
##  - attr(*, "class")= chr "classIntervals"
breaks_qt$brks
## [1]    0.00000   53.03633  134.37850  260.75619  470.58824 1343.87352
spplot(philly, "HOMIC_R", col.regions=pal, at = breaks_qt$brks,  main = "Philadelphia homicide rate per 100,000")

# add a very small value to the top breakpoint, and subtract from the bottom for symmetry 
br <- breaks_qt$brks 
offs <- 0.0000001 
br[1] <- br[1] - offs 
br[length(br)] <- br[length(br)] + offs 

# plot
spplot(philly, "HOMIC_R", col.regions=pal, at=br,  main = "Philadelphia homicide rate per 100,000")

# correcting legend
philly$HOMIC_R_bracket <- cut(philly$HOMIC_R, br)
head(philly$HOMIC_R_bracket)
## [1] (261,471]   (-1e-07,53] (53,134]    (-1e-07,53] (-1e-07,53] (-1e-07,53]
## Levels: (-1e-07,53] (53,134] (134,261] (261,471] (471,1.34e+03]
class(philly$HOMIC_R_bracket)
## [1] "factor"
spplot(philly, "HOMIC_R_bracket", col.regions=pal, main = "Philadelphia homicide rate per 100,000")

# install.packages('GISTools')
library(GISTools)                               # load library
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
# View(philly)
choropleth(philly, philly$HOMIC_R)              # plot the polygons
shd <-  auto.shading(philly$HOMIC_R)            # we need that for the legend coloring
choro.legend(                                   # plot the legend
  bbox(philly)["x","max"] - 5000,               # x coordinate of top left corner
  bbox(philly)["y","min"] + 15000,              # y coordinate of top left corner
  shd                                           # color scheme
  )                               
title("Philadelphia homicide rate per 100,000") # plot the title.

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
# getwd
philly_sf <-  st_read("/Users/josephhaymaker/Desktop/STAT571_Denver_Animal_Protection/Philly3/Philly3.shp")
## Reading layer `Philly3' from data source `/Users/josephhaymaker/Desktop/STAT571_Denver_Animal_Protection/Philly3/Philly3.shp' using driver `ESRI Shapefile'
## Simple feature collection with 381 features and 12 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 1739498 ymin: 457288.9 xmax: 1764018 ymax: 490539.6
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
plot(philly_sf)
## Warning: plotting the first 10 out of 12 attributes; use max.plot = 12 to
## plot all
## Warning in classInt::classIntervals(values, nbreaks, breaks): var has
## missing values, omitted in finding classes

plot(philly_sf$HOMIC_R) # this is a numeric vector

plot(philly_sf["HOMIC_R"])

hr_cuts <-  cut(philly_sf$HOMIC_R, br)
plot(philly_sf["HOMIC_R"], main = "Philadelphia homicide rate per 100,000", col = pal[as.numeric(hr_cuts)])
legend(1760000, 471000, legend = paste("<", round(br[-1])), fill = pal)

# GGPLOT
#So if we wanted to map polygons, like census tract boundaries, we would use longitude and latitude of their vertices as our x and y values and  geom_polygon() as our geometry.

library(broom)
philly_df <- tidy(philly)
## Regions defined for each Polygons
dim(philly_df) #31410     7
## [1] 31410     7
# head(philly_df, 30)
# View(philly_df)
# head(philly_df)
str(philly, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  381 obs. of  13 variables:
##   ..@ polygons   :List of 381
##   ..@ plotOrder  : int [1:381] 378 377 11 182 25 16 372 3 94 360 ...
##   ..@ bbox       : num [1:2, 1:2] 1739498 457289 1764018 490540
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
philly$polyID <- sapply(slot(philly, "polygons"), function(x) slot(x, "ID"))
philly_df <- merge(philly_df, philly, by.x = "id", by.y="polyID")
head(philly_df)
##   id    long      lat order  hole piece group OBJECTID     AREA PERIMETER
## 1  0 1758028 490539.6     1 FALSE     1   0.1        1 37863520  32091.95
## 2  0 1758067 490520.2     2 FALSE     1   0.1        1 37863520  32091.95
## 3  0 1758134 490486.9     3 FALSE     1   0.1        1 37863520  32091.95
## 4  0 1758294 490407.7     4 FALSE     1   0.1        1 37863520  32091.95
## 5  0 1758369 490370.2     5 FALSE     1   0.1        1 37863520  32091.95
## 6  0 1758451 490329.6     6 FALSE     1   0.1        1 37863520  32091.95
##   TRACTID FIPSSTCO TRT2000       STFID N_HOMIC TOTALPO mdHHnc_  HOMIC_R
## 1     365    42101  036500 42101036500       3     703      NA 426.7425
## 2     365    42101  036500 42101036500       3     703      NA 426.7425
## 3     365    42101  036500 42101036500       3     703      NA 426.7425
## 4     365    42101  036500 42101036500       3     703      NA 426.7425
## 5     365    42101  036500 42101036500       3     703      NA 426.7425
## 6     365    42101  036500 42101036500       3     703      NA 426.7425
##   PCT_COL HOMIC_R_bracket
## 1 1.78117       (261,471]
## 2 1.78117       (261,471]
## 3 1.78117       (261,471]
## 4 1.78117       (261,471]
## 5 1.78117       (261,471]
## 6 1.78117       (261,471]
View(philly_df)

library(ggplot2)

ggplot() +                                               # initialize ggplot object
  geom_polygon(                                          # make a polygon
    data = philly_df,                                    # data frame
    aes(x = long, y = lat, group = group,                # coordinates, and group them by polygons
        fill = cut_number(HOMIC_R, 5))) +                # variable to use for filling
  scale_fill_brewer("Homicide Rate", palette = "OrRd") + # fill with brewer colors 
  ggtitle("Philadelphia homicide rate per 100,000") +    # add title
  theme(line = element_blank(),                          # remove axis lines .. 
        axis.text=element_blank(),                       # .. tickmarks..
        axis.title=element_blank(),                      # .. axis labels..
        panel.background = element_blank()) +            # .. background gridlines
  coord_equal()                                          # both axes the same scale

# ggplot(philly_sf) + geom_sf(aes(fill=HOMIC_R))
library(tmap)
names(philly)
##  [1] "OBJECTID"        "AREA"            "PERIMETER"      
##  [4] "TRACTID"         "FIPSSTCO"        "TRT2000"        
##  [7] "STFID"           "N_HOMIC"         "TOTALPO"        
## [10] "mdHHnc_"         "HOMIC_R"         "PCT_COL"        
## [13] "HOMIC_R_bracket" "polyID"
tm_shape(philly) +
  tm_polygons("HOMIC_R", style="quantile", title="Philadelphia \nhomicide rate \nper 100,000") +
  tm_shape(philly) + tm_polygons("TOTALPO", style="quantile") +


tmap_mode("view")
## tmap mode set to interactive viewing
last_map()
## tmap mode set to interactive viewing
vignette("tmap-nutshell")
## starting httpd help server ...
##  done
#  ============================================
tm_shape(denver4) +
  tm_polygons("NBHD_NAME", style="quantile", title="Denver county statistical neighborhoods")
# tmap_mode("view")
# last_map()

tm_shape(denver5) +
  tm_polygons("DISTRICT_I", style="quantile", title="Philadelphia \nhomicide rate \nper 100,000")